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We construct a nonlocal density functional approximation with full exact exchange, while pre- 
serving the constraint-satisfaction approach and justified error cancellations of simpler semilocal 
functionals. This is achieved by interpolating between different approximations suitable for two 
extreme regions of the electron density. In a "normal" region, the exact exchange-correlation hole 
density around an electron is semilocal because its spatial range is reduced by correlation and be- 
cause it integrates over a narrow range to —1. These regions are well described by popular semilocal 
approximations (many of which have been constructed nonempirically) , because of proper accu- 
racy for a slowly-varying density or because of error cancellation between exchange and correlation. 
"Abnormal" regions, where nonlocality is unveiled, include those in which exchange can dominate 
correlation (one-electron, nonuniform high-density, and rapidly- varying limits), and those open sub- 
systems of fluctuating electron number over which the exact exchange-correlation hole integrates 
to a value greater than —1. Regions between these extremes are described by a hybrid functional 
mixing exact and semilocal exchange energy densities locally, i.e., with a mixing fraction that is a 
function of position r and a functional of the density. Because our mixing fraction tends to 1 in the 
high-density limit, we employ full exact exchange according to the rigorous definition of the exchange 
component of any exchange-correlation energy functional. Use of full exact exchange permits the 
satisfaction of many exact constraints, but the nonlocality of exchange also requires balanced non- 
locality of correlation. We find that this nonlocality can demand at least five empirical parameters, 
corresponding roughly to the four kinds of abnormal regions. Our local hybrid functional is perhaps 
the first accurate fourth-rung density functional or hyper-generalized gradient approximation, with 
full exact exchange, that is size-consistent in the way that simpler functionals are. It satisfies other 
known exact constraints, including exactness for all one-electron densities, and provides an excellent 
fit to the 223 molecular enthalpies of formation of the G3/99 set and the 42 reaction barrier heights 
of the BH42/03 set, improving both (but especially the latter) over most semilocal functionals and 
global hybrids. Exact constraints, physical insights, and paradigm examples hopefully suppress 
"overfitting" . 

PACS numbers: 31.15.ej, 31.10.+Z, 71. 15. Mb 



I. INTRODUCTION 

Kohn-Sham density functional theory [l|, 0] is a com- 
putationally efficient and often usefully accurate ap- 
proach to electronic structure calculation in condensed 
matter physics and quantum chemistry. In this the- 
ory, the ground-state electron density and total energy 
E are found by self-consistent solution of a one-electron 
Schrodinger equation. Only the exchange-correlation 
(xc) energy as a functional of the density needs to be 
approximated in practice. This energy can always be 
written as 

£7 xc [n T ,nj] = / d 3 rn(r)e KC {[n h n l ];r). (1) 



In this equation, the energy density is the product of 
the electron density n(r) = rij(r) + rij(r) (the sum of 
the spin densities) and the position-dependent exchange- 
correlation energy per electron e xc ([n-f , nj; r). 

Local or semilocal approximations for e xc ([ri|, nj; r), 
as defined below, are computationally fast and concep- 
tually simple, and can be constructed without empiri- 
cism. They work best for sp-block atoms and their 
molecules and solids around equilibrium. In real sys- 
tems, full nonlocality is most needed to correct the so- 
called self-interaction errors [3J, which can be most se- 
vere in certain stretched-bond situations that arise in 
transition states for chemical reactions and in the dis- 
sociation limit (and also effectively in the d- and /-block 
systems) 0, H, i, 0, H, i, ES H3 ■ In this paper, we follow 
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a conservative approach to full nonlocality: We identify 
those "normal" regions of space in which semilocality is 
justified, and those "abnormal" ones in which it is not, 
and then introduce full nonlocality only to the extent that 
a region is abnormal. The density functional that we con- 
struct here shows promising early results, but, whether 
or not that promise is ultimately fulfilled, the underlying 
insights could have lasting value. 

In popular approximations, e xc ([n-|-, nj; r) is a func- 
tion of ingredients at position r, and different selections 
of in gred ients define different rungs of a "Jacob's lad- 
der" [l2| of approximations. Accuracy tends to increase 
up the ladder, while simplicity increases downwards. The 
first three rungs use semilocal ingredients constructed 
from the density or orbitals in an infinitesimal neigh- 
borhood of r, and have nonempirical [l3| (as well as em- 
pirical) constructions. The first rung is the local spin 
density approximation (LSDA) [l|, 0, [TH, which uses 
as ingredients n^(r) and 7ij.( r ) and is exact only for a 
uniform or very slowly- varying electron gas. The second 
rung is the generalized gradient approximation (GGA), 
which adds the gradients Vn-|-(r) and Vn^(r), as in the 
Perdew-Burke-Ernzerhof (PBE) [jj] and PBEsol [l7j] 
nonempirical constructions. The third rung is the meta- 
GGA, which adds the positive orbital kinetic energy den- 
sities Tf(r) and Tl(r), as in the Tao-Perdew-Staroverov- 
Scuseria (TPSS) [H, [13, H3] nonempirical construction. 
Alternatively, the Laplacians V 2 7ij(r) and V 2 n^(r) of the 
spin densities could be used [2l|. Self-consistent calcula- 
tions with any of the semilocal functionals arc relatively 
easy and efficient, especially if the demand for a single 
multiplicative effective potential is dropped at the meta- 
GGA level (as in Refs. Ill, [H HJ . 

The fourth rung or hyper-GGA 0, H H, [H H [H, 
US US US H3] i with which we will be concerned here, adds 
a fully nonlocal (i.e., not semilocal) ingredient, the exact- 
exchange energy per electron e x x . Unlike total energies E, 
energy densities ne xc are non-unique or gauge-dependent. 
In the conventional gauge, 

ef(ron v )(r) = i/dV^^, (2) 

where 

, n ST^ l^r(r, r')| 2 
Mr,r) = -^ n(r) , (3) 

is the exact-exchange hole density around an electron at 
r, and 

P<7 {v,v')=Y,MU^W) (4) 

i 

is the Kohn-Sham one-particle density matrix for spin a. 
The fi a are occupation numbers, typically or 1 except 
for values in between that can arise when a Kohn-Sham 
orbital ^ ?;[T (r) is shared with another system [3, 0, [3, H, 
0, H, H, GjJ • The exact-exchange hole density satisfies 



the sum rule 

= -1 + A x (r), (0<A X <1) (5) 

and the right-hand side of Eq. (|5|) reduces to —1 when all 
occupation numbers are or 1. A rigorous proof of the 
exchange-only equations ©-([5]) for an open system of 
fluctuating electron number is given in Ref. ITol . There is 
also a fifth rung (e.g., Ref.l3ll) which adds all the occupied 
and unoccupied orbitals as ingredients. While the fourth- 
and fifth-rung functionals can also be implemented fully 
self-consistently, tests and applications using, say, meta- 
GGA orbitals (as in the present work) are easier and more 
efficient and often suffice. 

As we ascend the ladder, the additional ingredients 
can be used to satisfy additional exact constraints on 
E KC [ni ,n J (nonempirical approach [13, El), or to fit 
data better (empirical approach), or both. As we will 
see here, many additional constraints can be satisfied on 
the fourth r ung, including exactness for all one-electron 
densities (3, ll2Ll23j and full exact exchange [13], but for 
the first time some empiricism becomes unavoidable (al- 
though it can be avoided again on the fifth rung |3ll|). 

There are good reasons for these particular choices of 
ingredients. The local density suffices for a uniform elec- 
tron gas. The density gradients contribute to second- 
order gradient expansions [33, [S3, Hij] for densities that 
vary slowly over space, and can also satisfy other con- 
straints [ill. The positive orbital kinetic energy density 
is relevant to exchange [H, HH and in particular to its 
fourth-order gradient expansion .18, 3G] , and can be used 
to zero out the correlation energy density in one-electron 
regions [H, [37| • The Kohn-Sham orbitals needed to con- 
struct the higher-rung ingredients are themselves fully 
nonlocal but implicit functionals of the density. (Strictly, 
the orbitals are density functionals when the optimum ef- 
fective potential is required to be a function of position r. 
In practice, ground-state densities and energies are not 
much affected when this requirement is dropped, as it 
often is for computational convenience in self-consistent 
implementations of meta-GGA's.) 

The ladder classification is not intended to be exclu- 
sive. Range-separated hybrids [3, H3, H3] that use the 
exact- exchange hole density of Eq. (J3J) fall slightly above 
the fourth rung but well below the fifth. Fully nonlocal 
"two-point" approximations or six-dimensional integrals 
that are explicit functionals of the density arc also pos- 
sible, and can provide a long-range van der Waals in- 
teraction [iflj ] that is still missing from our local hybrid 
functional, but implicit density functionals with a non- 
local dependence upon the orbitals are probably needed 
to describe the situation in which electrons are shared 
between open subsystems 0, i, @, 0, i, H, [3 • 

There are two principal reasons why density function- 
als should be constructed to satisfy known exact con- 
straints. The first is idealistic: Our approximations 
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should not needlessly violate what we know to be true. 
The second is practical: Satisfaction of many constraints 
helps the approximation to work over a wider range of 
densities. For example, functionals that are constructed 
only by fitting to molecular data are less accurate for 
solids than functionals that satisfy solid-state-like con- 
straints such as the uniform gas limit [4l| , which by itself 
fully defines the first or LSDA rung of the ladder. 

At the second or GGA rung, the constraints are al- 
ready less constraining. There are several different sets 
of constraints that can be imposed on this rung, and 
some are (for this limited form) incompatible with oth- 
ers [IH, [T3|. Even for a chosen set of constraints, there 
may be quite different ways in which they can be satis- 
fied . Thus, the nonempirical GGAs are in practice 
guided by some additional physical postulate. For the 
PBE GGA, this guiding postulate is the sharp (atomic- 
like) cutoff of the spurious long-range part of the gradient 
expansion for the exchange-correlation hole [43[ . For the 
PBEsol GGA for solids, it is restoration of the gradient 
expansion for the exchange energy over a wide range of 
slowly- or moderately- varying densities [13] • 

We close this introduction by summarizing some of 
the successes and failures of the TPSS meta-GGA, the 
nonempirical functional on the highest semilocal rung 
of the ladder. For the standard enthalpies of forma- 
tion of the 223 molecules of the G3/99 test set 0] (in- 
cluding COF2), the mean absolute error (MAE) com- 
puted self-consistently using the fully uncontracted 6- 
3f 1++G(3df,3pd) basis set is only 6.5 kcal/mol when 
TPSS exchange is combined with TPSS correlation. The 
MAE increases to 30.3 kcal/mol when full exact ex- 
change is combined with TPSS correlation, and to 211.0 
kcal/mol for exact exchange and no correlation (i.e., 
Hartree-Fock) . Thus, atoms and molecules at equilib- 
rium, with integer electron numbers, are predominantly 
normal systems in which the errors of semilocal exchange 
and semilocal correlation tend to cancel. This does not 
mean that semilocal functionals always predict realistic 
binding energy curves. In fact, they do so only when the 
fragments have integer electron numbers. Asymmetric 
molecules, when described by semilocal functionals, of- 
ten dissociate to fragments of spurious fractional charge, 
for which the energy is seriously too low (Ref. @ and ref- 
erences therein). Even when the fractional charge is real, 
as in the dissociation of Xj, the self-interaction error 
of semilocal functionals makes the energy seriously too 
low (Ref. [D, and references therein) at the dissociation 
limit X +1 / 2 • • ■ X +1 / 2 , where it ought to be equal to that 
of X° • • • X +1 . The combination of stretched bonds and 
noninteger average (hence fluctuating) electron numbers 
on the fragments makes for a highly abnormal valence re- 
gion. A milder version of the same abnormality is found 
for the moderately stretched bonds in the transition state 
of a chemical reaction, where semilocal functionals pre- 
dict reaction barriers that are too low or even negative. 



II. LOCAL HYBRID FORM OF THE 
HYPER-GGA 

The earliest hyper-GGAs were the global hybrid func- 
tionals [ijj Hll |4y, [47| , for which the simplest form is 

E* = aE? + (l-a)E* + E?. (6) 

Here E^ and E^ 1 are semilocal exchange and correlation 
functionals, E° x is the exact-exchange energy, and a is 
a universal, position-independent empirical mixing frac- 
tion between and 1. While semilocal approximations 
(a = 0) typically overestimate the atomization energies of 
molecules [ill |48[ and underestimate the energy barriers 
to chemical reactions (49|, the use of exact or Hartree- 
Fock exchange (a = 1) makes opposite errors. Fitting 
Eq. ([6]) to atomization energies often gives a rs 0.2, while 
the reaction barrier heights typically require a 0.5. 
For some systems and properties, even a « 0.2 is too 
large [5p| . The smallness of a indicates that semilocal 
exchange is more compatible with semilocal correlation 
than is the exact, fully nonlocal exchange. The global hy- 
brid functionals are perhaps the most popular functionals 
in modern quantum chemistry. (Screened hybrid func- 
tionals [5 if are also gaining popularity for solid-state cal- 
culations [52j, |53j). They often improve upon the semilo- 
cal functionals, at the sometimes-acceptable additional 
cost of a Hartree-Fock calculation, and the improvement 
can be pushed further by the addition of further empiri- 
cal parameters 47]. However, they satisfy no exact con- 
straint beyond those satisfied by the underlying semilocal 
E*. 

A natural generalization to Eq. © is the local hybrid 
functional 

4 h c (r) = o(r)^(r) + [l-a(r)]^(r)+^(r) 

= ef(r) + [1 - a(r)][4'(r) - ef (r)] + ef(r)(7) 

where < a(r) < 1. Eq. |7]) was first proposed in Ref. |22L 
without a form for a(r). Forms were proposed in Ref. 1 121 
and Ref. l23j (where the term "local hybrid" was coined), 
but those forms did not achieve high accuracy for equi- 
librium properties of molecules. 

The choice = eJ c PSS and a(r) = 1 in Eq. ([7]) satisfies 
nearly all exact constraints that a hyper-GGA can satisfy, 
but yields very poor atomization energies of molecules 
(as shown in Ref. [12, for example) because it misses the 
delicate and helpful error cancellation between semilocal 
exchange and semilocal correlation that typically occurs 
(because the exact xc-hole is deeper and more short- 
ranged than the exact x-hole) in "normal" regions of 
space (as defined in section UTTj) . From this fact, it is 
already clear that the known exact constraints do not 
tell us how much full nonlocality is needed in correlation 
for compatibility with exact exchange. The imbalance 
between nonlocal approximate exchange and semilocal 
approximate correlation may also be responsible for the 
inaccuracies of the Perdew-Zunger self-interaction cor- 
rection 0] and of the local hybrid functional of Ref. [23l . 
when applied to atomization energies. 
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As envisioned in Ref. [Tl and discussed below in sec- 
tion |IV3 a hyper-GGA has full exact exchange if it sat- 
isfies the exact constraint [H3| (for any system with a 
nondegenerate Kohn-Sham noninteracting ground state) 



lim 



1 



under uniform density scaling 

n(r) — » n\(r) = A 3 n(Ar) 



(8) 



(9) 



to the high-density limit A — > oo. (Note that n(r) and 
n^(r) have the same electron number). Such a func- 
tional will automatically satisfy all exact constraints on 
exchange alone. Of several developed hyper-GGAs 0, 
HI, HIM HE IH [Hj], the only one that seems to have 
this property is that of Refs. [IE H3, and HH, for which 
the correlation energy is not exactly size-consistent. Our 
size-consistent local hybrid has full exact exchange be- 
cause its mixing fraction a(r) tends rapidly enough to 
1 in the high-density limit, while e° x ~ A, ~ A, and 
e^ 1 ~ A in this limit. Formally, the exchange part of any 



is m 



„ r i ,. E xc [n\] 
rj x [n\ = lira 

A — >oc A 



(10) 



Equation (10) is the only non-arbitrary way to define 
the exchange component of a density functional for the 
exchange-correlation energy. 

Thus, the three terms in the last line of Eq. ([7]) rep- 
resent respectively exact exchange, the fully nonlocal 
part of correlation including "left-right" [56| static or 
near-degeneracy correlation important in molecules, and 
the semilocal dynamic correlation. Instead of thinking 
of a(r) as the fraction of exact exchange mixed locally 
with semilocal exchange, it is more correct to think of 
[1 — a(r)] as the fraction of the nonlocal energy difference 
[e^(r) — £x X ( r )] added to exact exchange plus semilocal 
correlation. The quantity [e^ ( r ) — e x X ( r )] ^ s ^ ne nonlocal- 
ity error (and arguably often the self-interaction error) of 
the semilocal e^ 1 - In a normal region (a « 0), it provides 
an estimate of the nonlocal part of the correlation energy 
(such as the static correlation). In an abnormal region, 
the semilocal functionals often overestimate [1(3] the non- 
local correction to correlation, and their estimate must 
be scaled back. Our a(r) is mnemonic for and measures 
the "abnormality" of a region of space. 

We have recently argued [l(j that Eq. (O could be 
useful not only for molecules and solids near equilibrium, 
but even for the more challenging problem of stretched 
bonds between open subsystems with noninteger average 
electron number. In such subsystems, in the limit of infi- 
nite bond stretch, the total energy variation with electron 
number between adjacent integers is concave downward 
for exact exchange [IJJ] , concave upward for semilocal ex- 
change or exchange-correlation [5(, and linear in an ex- 
act description justifying the local- hybrid mixing of 
Eq. (O . To the extent that a functional mimics the exact 



linearity, it is said to be "many-electron self-interaction- 
free" [1, 0, [|| . We have also discussed [l(| in a general 
way how a(r) should be much less than 1 in normal re- 
gions where ej c can be accurate, but should approach 1 
in many abnormal regions where e^ c cannot be accurate. 
We shall construct a(r) in sections HTT1 and [TVl This con- 
struction will necessarily be complicated, because we aim 
to satisfy as many exact constraints as possible. Our con- 
struction will explain why the more complicated hyper- 
GGAs [13, HE US H3 (including ours) require typically 
four or five empirical parameters to balance the full non- 
locality of correlation with that of exchange. (The hyper- 
GGAs of Refs. H3, [H, [H, H3 have additional empirical pa- 
rameters in their semilocal parts, while our hyper-GGA 
does not.) 

To satisfy as many constraints as possible, we will take 



4'(r)= ex TPSS (r), ^(r)= £ J PSS (r). 



(11) 



The TPSS meta-GGA [3 EH El is constructed nonem- 
pirically to satisfy those exact constraints that a meta- 
GGA can, including the constraint e c (r) = in one- 
electron regions. 

A complication arises because of the non-uniqueness 
of the exchange energy density: many different energy 
densities can integrate to the same energy. In Eq. 0, 
we want eTUr) to be as similar to e^fc) as it can be, so 
we write [571 



n(r)4 x (r) = n(r)e^ conv )(r) + V • [f/(r)W(r)], (12) 

where e ® x ( conv ) j s the conventional exact-exchange energy 
per electron from Eq. ([2]) . The functions U (r) and V(r) 
are constructed from hyper-GGA ingredients in Ref. [57|. 
The divergence of the vector field U W in Eq. (fT2")) inte- 
grates over r to zero. The divergence term is largest in 
regions where the density is dominated by a single orbital 
shape, such as one-electron and tail regions; omitting it 
increases the errors of our local hybrid functional quan- 
titatively but not qualitatively. The global hybrids of 
Eq. are of course invariant under any choice of gauge. 

We have argued that the empirical and nonempiri- 
cal approaches to density functional approximation must 
flow together at the hyper-GGA level, as they need not on 
the three lower rungs of the ladder. On this fourth rung, 
the world is almost turned upside down: Up to a point, 
building more physics into a hyper-GGA demands more 
and not fewer empirical parameters. (A similar situation 
arises in the construction of explicit density functionals 
for the orbital kinetic energy |21j.) 

Simple one-parameter local hybrid functionals, satis- 
fying few or no exact constraints beyond those of the 
lower-rung functionals, have been constructed by Kaupp 
and collaborators [H, H^]- More complicated function- 
als with four or five empirical parameters have been de- 
veloped using all three standard physical approaches to 
density functional approximation: 

(a) Modeling [35|, [43| the exchange-correlation hole 
around an electron, including the static-correlation part 
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that arises in stretched-bond H2, leads to Becke's hyper- 
GGA [2H [H| . Because Becke starts from a system that 
has an exact degeneracy of its Kohn-Sham noninteract- 
ing ground state, he is able to avoid some of the symme- 
try breaking that persists in our and other hyper-GGAs, 
but at the cost of violating the exact constraint of Eq. ([8|) 
for all ground states without such an unusual degener- 
acy. For any system in which Becke's hyper-GGA and 
ours predict a non-zero static correlation energy, the ra- 
tio of static correlation to exchange energies under uni- 
form density scaling to the high-density limit will remain 
unchanged for Becke's hyper-GGA while it will tend to 
zero for ours; the latter behavior is the correct one in the 
normal case where the static correlation in the unsealed 
system arises from a near- but not exact degeneracy. 

(b) Modeling [II, [H, HI] the adiabatic connection [59| 
between the Kohn-Sham Slater determinant and the true 
interacting wave function leads to the hyper-GGA of 
Yang and collaborators [26|, [27|, but with a loss of size 
consistency (in the weak sense in which simpler func- 
tional are size-consistent). Without size consistency, the 
energy of a subsystem (e.g., an atom or molecule) can de- 
pend upon the presence of another (e.g., a metal surface) 
even in the limit of infinite separation and even without 
charge transfer between the subsystems. Size consistency 
can be preserved by using energy densities, as here and 
in Refs. [U, HH. It is lost when total energies are used 
nonlinearly. 

(c) Satisfaction [l6|, [l^, [6(| of exact constraints on 
E^ c [n^,ni] leads to our local hybrid. 

Hyper-GGAs have also been constructed by fitting a 
large number of empirical parameters to chemical data 
sets (e.g., Refs . l47l and l6lh . while satisfying relatively few 
exact constraints. In particular, the M06-HF hyper-GGA 
of Zhao and Truhlar [61| is claimed to have "full Hartree- 
Fock exchange" , but like Becke's hyper-GGA [H, [25| it 
contains in addition other terms (semi-local ones in this 
case) which improperly scale like exchange in the high- 
density limit. Plots [47| of the enhancement factors of the 
semi-local parts of the parent functional M05-2X display 
wiggles and other anomalous behaviors suggesting more 
parameters and fewer constraints than are needed. 

Our hyper-GGA work here is distinguished from that 
of others in two principal ways: (i) We aim to under- 
stand fully and build in the error cancellation [Iol | be- 
tween semilocal exchange and semilocal correlation that 
is responsible for much of the success of simpler func- 
tionals. (ii) We aim to satisfy nearly all the known exact 
constraints that a hyper-GGA can, in order to test and 
apply this exact knowledge and to make our functional 
useful even in systems and situations very different from 
those where we have fitted or tested it. One example 
of such a system is jellium uniaxially compressed toward 
the true two-dimensional limit [62| . For this challenging 
problem of dimensional cross-over, our local hybrid func- 
tional is significantly better than LSDA, PBE GGA, or 
TPSS meta-GGA, although not as good as the fifth-rung 
functional described in Ref. |3l|. 



III. WHY THE EXACT 
EXCHANGE-CORRELATION HOLE IS 
SEMILOCAL IN NORMAL REGIONS 

The exact exchange-correlation energy can be writ- 
ten [H, [63[ in the form of Eq. (TT]), where in the con- 
ventional gauge 

£ e,(co„v )(r) = 1 J d3r , J\ a n%«(r,r>) 
If 00 

= - duinun^r, u). (13) 
2 Jq 

Here n x x (r, u) is the spherical average over the direction 
of u = r' — r of the average over coupling constant a of 
n™' a (r,r'), where the density is held fixed as a varies 
from to 1 in the Coulomb interaction a/\r' — r|, and 
n xc'" =0 ( r ; r ') = n x X ( r J r ')- The exact xc-hole obeys the 
sum rule 

/>OC 

/ du±-KU 2 n™{r,u) =-l + A xc (r), (0 < A xc < 1), 
Jo 

where A xc vanishes if the integrated electron number 
in the system does not fluctuate [f| [l(|. In the ab- 
sence of exact degeneracy of the Kohn-Sham noninter- 
acting system, the a = limit of — 1 + A xc is shown 
on the extreme right-hand side of Eq. ([5]). For the 
slowly-varying densities on which the semilocal approx- 
imations are based, A xc = A x = 0. For the LSDA, 
n x f DA (r, u) = n x " lf (rif (r), n^(r); u); for the other semilo- 
cal approximations, the hole can be either modeled from 
the start [H, HH or reverse-engineered [H, [H[ . 

Gunnarsson and Lundqvist [631 ] gave an argument for 
the success of LSDA and hence of other semilocal approx- 
imations: These approximations satisfy the correct sum 
rule (when A xc = 0), which is by Eq. (fl4|) the second 
moment of the hole density. Then the energy per elec- 
tron, which is by Eq. (fT3|) the first moment of the hole 
density, should not be too wrong. This argument was ex- 
panded by Burke, Perdew and Ernzerhof [661 ] who argued 
that the small-u behavior (through order u 2 ) of « x x (r, u) 
is determined exactly by semilocal information, and the 
small- u behavior (through order \u\) of n™(r, u) is often 
determined approximately by local information. They 
also argued that the system average arising in Eq. ^ 
unweights regions of space in which the small-u behavior 
of n. x x (r, u) is not so well described by semilocal informa- 
tion. This system average also eliminates the differences 
between different energy-density gauges. 

This argument is easily extended to explain why the 
semilocal approximations usually work better for ex- 
change and correlation together than for either sepa- 
rately: Correlation makes n x x (r,it) deeper at small it, 
and more short-ranged, than n x x (r,w). A deep, narrow 
xc-hole density, integrating to —1 and having a semilocal 
small-w behavior, is semilocal. 

We define a normal region of space [10( as one in which 
n x x (r, u) is modeled reasonably well by a semilocal ap- 
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proximation, at least after the system average over that 
region. The electron density in a normal region is either 
reasonably slowly-varying, or reasonably low and many- 
electron-like, or both. Only at sufficiently low and many- 
electron-like density is correlation comparable in strength 
to exchange, making possible a useful error cancellation 
between them. In a normal region, we want a(r) « 
in Eq. to take advantage of the proper accuracy of 
semilocal approximations (for a slowly-varying density) 
or of the error cancellation between semilocal exchange 
and semilocal correlation. 

Semilocal information (n, Vn, r) is enough to tell 
us whether a density is slowly-varying, low, or many- 
electron-like. We shall take advantage of this in sec- 
tion IIV Al The fully nonlocal information in e x x of 
Eqs. @ and (TT2]) can tell us how big A x is. Exact 
or near degeneracies can make A xc small or zero even 
when A x is not small or zero (as in infinitely-stretched 
spin-unpolarizcd H 2 ), but in those cases spin-symmetry 
breaking can still rescue the semilocal a ppr oximations 
by restoring the correct small-it behavior [67}. We shall 
make use of this in section HV Bl 

In a normal region, we will take a(r) in Eq. ([7]) to be 
small or zero. Conversely, an abnormal region is one in 
which we have no clear reason to keep a(r) small, and 
can let a(r) approach 1. 

By the way, we can write the spherically-averaged 
exchange-correlation hole that yields our local hybrid 
functional of Eq. (UJ as 



..Hi 



(r,u) = n™(r,u) + [l-a(r)][n^(r,u) 



n x x (r, u)} 



TPSS 



(r, u). 



(15) 



The TPSS hole is known [65j. Note that the middle term 
on the right of Eq. (|T5)) is purely long range in u [and of 
infinite range when n° x (r,w) is], confirming our interpre- 
tation of this term as the nonlocal or static correlation. 



to identify an abnormal region according to conditions 
(i) or (ii), respectively. We will then combine these into 
a single a(r) in subsection IIV CI 



A. Abnormal regions where exchange can 
dominate: one-electron, rapidly-varying, and 
nonuniform high-density regions 

As discussed earlier, no error cancellation between ex- 
change and correlation can be expected in abnormal 
regions where the density is too one-electron-like, too 
rapidly- varying over space, or too high, so there is no rea- 
son not to use full exact exchange there. In regions where 
the density is slowly-varying over space, we might accu- 
rately use either full exact exchange (a = 1) or semilocal 
exchange (a = 0), but semilocal exchange is computa- 
tionally preferable because it avoids the need to integrate 
over the long tail of the exact exchange hole. 

Consider the interesting density- and position- 
dependent variable 



where 



GL2TPSS 



GL2TPSS 


,LSD ■ 



TPSS 



(16) 



([n TA) n u ];A 1 r) 



= lim e' c 

X — >oc 

is the Gorling-Levy [68j second-order or high-density 
limit of the TPSS correlation energy per electron at po- 
sition r. The latter tends , IT5 , to a negative finite 
limiting function of Ar as A — > oo under the uniform den- 
sity scaling of Eq. ©. The A -1 factor in Eq. (HJj then 
restores the length scale of the original density n(r). An 
explicit formula for e GL2TPSS j g gj ven j n Appendix [A"l 
Moreover, 



IV. MIXING EXACT EXCHANGE IN 
ABNORMAL REGIONS WHERE THE HOLE IS 
FULLY NONLOCAL 

Our hyper-GGA is the local hybrid functional of 
Eqs. 0, ©, and ([TT]). The motivation for this form 
was presented in section [TlJ and the idea behind the lo- 
cal exact-exchange mixing fraction a(r) was explained in 
section IIIII and Ref. 10: We want to satisfy essentially 
all possible exact constraints on E xc [n^, nj, including 
Eq. {5} which guarantees that our hyper-GGA has full 
exact exchange, while preserving the proper accuracy or 
the error cancellation between semilocal exchange and 
semilocal correlation that occurs in normal regions 

Thus we want a(r) ~ in a normal region, and 
a(r) ~ 1 in strongly abnormal regions. An abnormal 
region is one in which the density is (i) one-electron-likc, 
rapidly-varying over space, or nonuniform and high, or 
(ii) strongly fluctuating in electron number. We will de- 
fine an ai(r) and an a 2 (r) in subsections IIV Al and IIV Bl 



LSD 



([n h ni];r) = -- 



3 /3\ 1/3 (2n T ) 4 / 3 + (2n x ) 4 / 3 



2 («T + n l) 



(18) 

is the LSD exchange energy per electron at r. Note that 
u vanishes in the one-electron and rapidly-varying limits 
(because £ G L2TPSS does), and u also vanishes in the high- 
density limit (because l/e x SD does). 

Now we choose a mixing fraction ai(r) as a function 
of u which falls monotonically from 1 at u — to as 
u — > oo. Experience with global hybrid functionals (in 
which a is independent of r) suggests that cii might be a 
weak function of u over a wide range of u. Thus we try 



1 



ai = 



1 + Aln(l + J5u)' 



(19) 



where A and B are positive empirical parameters. This 
function is plotted in Fig. [TJ 

In Eq. ([TJ, the first term on the extreme right is the 
exact-exchange energy per electron, while the sum of the 
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FIG. 1: Mixing fraction ai(u) of Eq. l[l9]). The values of 
parameters A and B were fitted as explained in Section IVl 



middle and last terms is the local hybrid correlation en- 
ergy per electron ej, h . We will now verify that many of 
the exact constraints satisfied by TPSS correlation are 
also satisfied by our local hybrid correlation. 

For a one-electron density, properly ej pss = [H, Hi| , 
so £ GL2TPSS = o, u = 0, and oi = 1. Thus properly 
- 11 1 — 0. Unlike the semilocal functionals, our local hybrid 



is exact for all one-electron densities (H, H^, etc.) 
In the rapidly- varying limit, in which 



|Vn| 



2(37T 2 ) 1 /3 n 4/3 



(20) 



/TPSS 



o e m 

Jh 



so e 



GL2TPSS 



0, 



0, and 



di — > 1. Thus, e" 1 — > 0. In the exponential tail of an 
electron density, this is the correct behavior. Under one- 
dimensional density scaling to the true two-dimensional 
limit [63, [Z0, Ell > 4 h is correctly finite, although in this 
case the limiting function should properly be negative 
and not zero. 

In the nonuniform high-density limit [A — > oo under 
the uniform scaling of Eq. ©], in which 0, we 

find u — > and a\ — > 1 — ABu. Then 



ABe GL2TPSS 



.TPSS 



C LSDA 



OL2TPSS 



(21) 



The limit on the right of Eq. (|2~lj) is a function of Ar, but 
one that does not otherwise scale with A. This is properly 
a finite limit [H[ (if not always a properly negative one), 
implying Eq. (JSJ) and demonstrating that our local hybrid 
has full exact exchange. The first term on the right of 
Eq. (f2Tj) can be interpreted as the Gorling-Levy limit of 
the static correlation term in our local hybrid correlation 
of Eq. ([7]) . (Only when the static correlation energy arises 
from an exact degeneracy of the Kohn-Sham nonintcr- 
acting system can the corresponding limit of the exact 
static correlation energy not be finite [To|.) 

There are two regions of an atom or molecule in which 
ax — ► 1 is expected (and found by our local hybrid): 



the rapidly- varying exponential density tail and the high- 
density region of the deep core. For the rapidly-varying 
density of a quantum well whose thickness shrinks down 
to zero (the true two-dimensional limit), a\ — > 1 is also 
expected (and found); our local hybrid improves greatly 
but imperfectly [62j upon TPSS and other semilocal func- 
tionals in this limit. 

In the slowly- varying limit, in which s and other di- 
mensionless density derivatives tend to zero, we have 



TPSS 



C LSD 



K 



IS, 1690, eJf L2iPbb -> -oo like Ins and 
u — > oo, so a\ — > like l/ln(— Ins) and ej, h — > ej pss . 
Since the TPSS correlation energy recovers [lj| [IH, [69[ 
the second-order gradient expansion [33| in this limit, so 
does our local hybrid functional. 

Finally, we note that in the many-electron low-density 
limit [A — » under the uniform scaling of Eq. ©], 
l/ e x SD -» _0 °j so u — > oo, clx — ► 0, and — > eJ c PSS - 
There are reasons [72j to believe that TPSS is accurate 
(although not exact) for the exchange-correlation energy 
in this limit: As correlation becomes stronger in com- 
parison with exchange, more error cancellation between 
TPSS exchange and TPSS correlation is expected. 

We have not been able to prove, for all possible den- 
sities, that the local hybrid correlation energy is al- 
ways negative or that the Lieb-Oxford bound on its 
exchange-correlation energy is always satisfied. These 
constraints are easier to guarantee with semilocal func- 
tionals [ll|[l8:, 69] than with local hybrids. 

This ai makes our exchange-correlation functional ex- 
act in the one-electron and high-density limits, and more 
correct in the rapidly- varying limit. It should provide 
a — 1/r asymptote for the effective exchange-correlation 
potential around an atom or molecule in the limit r — ► oo, 
needed to bind small negative ions [3[ and also needed in 
time-dependent density functional theory when the 
electron density is driven far from the nuclei. But it can- 
not deal with the class of problems addressed in the next 
section. 



B. Abnormal regions where the hole density does 
not integrate to — 1: Open subsystems of fluctuating 
electron number 

Semilocal approximations to the density functional for 
the exchange-correlation energy assume that the exact 
exchange-correlation hole density around an electron in- 
tegrates to —1 on the scale of the local Fermi wavelength, 
as it does in an electron gas of slowly-varying density. 
Regions where this is not true, e.g., open subsystems 
of fluctuating electron number where A xc of Eq. (fT?)) is 
nonzero 0, [13] , are therefore abnormal. 

Consider the second interesting density- and position- 
dependent variable 



e TPSS ' 



(22) 



Where v is sufficiently less than 1, we can conclude [by 
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M 1/2 = ig/rf 

FIG. 2: Function f(v) of Eq. ([23}. The value of C was found FIG. 3: Function g{x) denned by Eq. J37). The values of 
as explained in Section [V] F is related to C by Eq. (|24|l . parameters 7? and i? were fitted as explained in Section IVl 



the argument around Eqs. (fTB")) and (fl"4"l) ] that A x of 
Eq. §5§ is positive. Because correlation suppresses den- 
sity fluctuation f73[ , we can expect A xc < A x . We seek a 
mixing fraction a 2 such that 02 ~ 1 when A xc « A x 3> 0, 
but ai « when A xc w 0. 
A first candidate for ai is 



r 1, «<c 

/(U)= l + cxp[l/(l-^-l/(,-C)T C<w<1 
I 0, w > 1 

where C < 1 is an empirical parameter. We have found 
that we can make f(v) flat around v — 1 or C, but not 
too steep at the midpoint v = (1 + C)/2, over a wide 
range of C < 1, by choosing 



F 



21n[(l-C)/2] 



> 0. 



(24) 



The function f(v) must be flat near v = 1 since values of 
i> slightly less than 1 can arise even in normal regions, ft 
interpolates smoothly between 1 (at the highly abnormal 
v < C) and (at the normal v w 1), reaching 1/2 at the 
midpoint v = (1 + C)/2, as illustrated in Fig. [2l If C is 
close to 1, there is a sharp step in f(v) which is smoothed 
when f(v) is multiplied by ej pss -e x x = (1 - v)e^ FSS . In 
practical calculations, the function f(v) is best evaluated 
not by Eq. (|23|) but as described in Appendix iBl 

If A xc were always equal to A x , we could take a 2 — 
f(v). But the inequality A xc < A x leads us to consider 
simple examples in which » < 1. This situation arises 
in certain stretched-bond diatomic molecules, with the 
bond length tending to infinity. 



(a) In stretched Hj [74j, a single electron is 



shared equally between two well-separated protons as 
H +1 /2 . . . H+ 1 / 2 . The electron fluctuates between the two 
proton centers, making A xc = A x > on each center. 
The Hartree-Fock or exact-exchange-only energy is cor- 
rectly equal to the sum of the Hartree-Fock energies for 



H and H+, while the TPSS meta-GGA energy is 0.102 
hartree (64.0 kcal/mol) too low compared to the sum of 
the TPSS energies for the integer-charge fragments. In 
this case, we want full exact exchange or 09 = 1. 

(b) In stretched neutral H2 [ToL [24 . 125I r75j , two elec- 
trons are shared equally between two well-separated pro- 
tons. In the exact singlet ground-state wave function, 
there is always exactly one electron on each proton cen- 
ter, because of static correlation arising from a Kohn- 
Sham degeneracy which becomes exact in the limit of in- 
finite bond length, and each H "atom" is spin-unpolarized 
as a result of fluctuation between "up" and "down" spin, 
making A xc = but A x > 0. The Hartree-Fock or exact- 
exchange-only energy is incorrectly 0.285 hartree (179 
kcal/mol) above the Hartree-Fock energy of two atoms. 
The corresponding TPSS error is reduced to 0.083 hartree 
(51.8 kcal/mol). In this case, we want no exact exchange 
or a 2 = 0. We can define an average v = 
where both F x x and E^ PSS are evaluated with the con- 
verged TPSS orbitals. For a predominantly normal sys- 
tem, v is close to 1 (e.g., 0.97 in Hj and 0.99 in sin- 
glet H2 at the respective equilibrium bond lengths), but 
v is significantly less than 1 in abnormal systems (0.62 
in infinitely-stretched H^ and 0.61 in infinitely-stretched 
singlet H 2 ). 

The hyper-GGA ingredient that seems to distinguish 
between these two cases is the relative spin polarization 



c 



(25) 



which equals ±1 for stretched H^ and for stretched H2. 
Thus we try 



«2 = 9 ( - ) f(v) 



where n = 3/47rr^ and, setting x = ( 2 /r 



9( x ) = 



Dx 



1 + Ex 



(26) 



(27) 



9 



in which D and E are positive empirical parameters, with 
D < E. Eq. (|27p is, like the exchange-correlation en- 
ergy itself, invariant under Q — > — C- It interpolates be- 
tween (at x = 0) and D/E (as x — » oo), as shown 
in Fig. [3] Our a 2 vanishes when ( — > (as in singlet 
H2), but also when r s — > 00 (i.e., in the low-density 
limit where correlation can become as strong as exchange, 
increasing the possibility of error cancellation). With 
this choice, it may be possible to get an improvement 
in much of time-independent chemistry and condensed 
matter physics from a 2 alone, although cii is still needed 
for the situations described at the end of section IIV Al 
The selection of £ as an ingredient of a 2 is not altogether 
satisfactory, since it is motivated by simple examples and 
not by general physical principles. 

Our simple examples will now be discussed further: (i) 
In stretched H 2 , we have a one-electron density which is 
already properly described by a\ = 1 from section HV Al 
(ii) In stretched H2, the standard spin-symmetry break- 
ing of semilocal functionals, which has a good physical 
basis [gl]], will localize a spin- up electron on one proton 
center, and a spin-down electron on the other, making 
\C\ = 1 almost everywhere and eliminating fluctuations 
(making A x = 0). 

Under the uniform density scaling of Eq. v(r) — > 
v(Xr) and £(r) — » C(Ar) an d ^s(Ar) — > A _1 r s (r). [Simi- 
larly, for s of Eq. (|20[) . s(r) — * s(Ar)]. In the high-density 
limit, our a 2 (like our a{) properly has an expansion in 
powers of r s (or A -1 ). 

We have invoked five empirical parameters: two (A 
and B) in a\ and three (C, D, and E) in a 2 . Next we 
will combine our two mixing coefficients into one. 

C. Combination rule for the mixing fractions 01 

and a,2 

We must combine the mixing fractions a\ and a 2 from 
sections llV Al and IIV Bl to get a single mixing fraction a 
for use in Eq. ([7]), with the following features: (i) < 
a < 1, as expected for any local or global hybrid, (ii) 
a = 1 when either a\ — 1 or 02 = 1, because either 
condition indicates a strongly abnormal region in which 
full exact exchange can be used, (iii) a = a 2 when a± = 
and a — a\ when 02 = 0, since in either case all the 
abnormality of a region is of either one type or the other. 
Condition (iii) also ensures that a = when a± = 02 = 0, 
indicating a strongly normal region in which semilocal 
approximations suffice. More generally, a > max(ai,a2), 
because a should increase with abnormality. 

A combination rule which satisfies all these expecta- 
tions is 

a = 1 — (1 — <zi)(l - a 2 ) = ai + a 2 - aia 2 . (28) 

Note that, in a nearly-normal region where a\ and a 2 
are small, we find a w a± + a 2 . Combining a± with a 2 
does not lose any of the exact constraints satisfied by a\ 
alone. The flatness of f(v) at v = 1 preserves the correct 
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FIG. 4: PSTS hyper-GGA mixing fractions a constructed for 
the Mg atom and isoelectronic Si 2+ and Ca 8+ ions using self- 
consistent TPSS orbitals obtained in Partridge's uncontracted 
basis sets [H, S El El : (20s,12p) for Mg; (20s,15p) for Si; 
(23s, I6p) for Ca. The parameters of a are given in Section IVl 
Here £ = 0, so 02 = and a = ai. The figure illustrates that 
our mixing fraction a tends smoothly to 1 in the low-density 
tail, and slowly to I in the core in the high-density (Z — » 00) 
limit. 



gradient expansion for correlation in the slowly-varying 
limit. The AB term of the high-density (A — > 00) limit 
of Eq. (|2~lj) is multiplied by the finite 1 — a2(Ar), and so 
remains finite. 

For the convenience of having a name, we will call 
the constructed local hybrid functional the PSTS hyper- 
GGA. The behavior of the PSTS hyper-GGA mixing 
fraction a and its components (with parameters fitted as 
described in Section [Vj) is illustrated by Figs.SHS] These 
figures are qualitatively reasonable, as discussed near the 
end of section IVII 



V. FITTING THE PARAMETERS TO 
STANDARD ENTHALPIES OF FORMATION 
AND BARRIER HEIGHTS 

We have determined the five empirical parameters 
(A, B, D, C, E) of the PSTS hyper-GGA mixing fraction 
o(r) by minimizing the quantity 

S w = i [MAE(G2/97) + MAE(BH42/03)] , (29) 

where MAE(G2/97) is the calculated mean absolute error 
of 148 standard enthalpies of formation of the molecules 
of the G2/97 test set 0, and MAE(BH42/03) is the 
mean absolute error of 42 forward and reverse barrier 
heights of the 21 gas- pha se hydrogen-transfer reactions of 
the BH42 /03 test set [ig| . The molecules of the G2/97 set 
are built up from hydrogen and/or first- and second- row 
atoms, and all but one of the reactions of the BH42/03 
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-6 -4 -2 2 4 6 8 
z (bohr) along the intemuclear axis 

FIG. 5: PSTS hyper-GGA mixing fraction a constructed with 
self-consistent TPSS/u-6-311++G(3df,3pd) orbitals along 
the intemuclear axis of a stretched F2 molecule (z — is 
the center of the bond). In all regions shown, 0,2 — and 
a — ai . The parameters of a are given in Section [V] Note 
that a has a cusp at z = which seems unphysical but is 
limited to a very small region of three-dimensional space. 
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FIG. 6: PSTS hyper-GGA mixing fractions a\ and a con- 
structed with self-consistent TPSS/u-6-311++G(3df,3pd) or- 
bitals along the intemuclear axis of a stretched NeJ molecular 
ion (z = is the center of the bond). The parameters of a 
are given in Section [V] In some regions, a 2 > and a > ox. 



test set involve atomic and molecular radicals (i.e., spin- 
unpaired systems) as reactants. The reference values of 
the standard enthalpies of formation are all experimental, 
while the 42 barriers are best estimates based on a combi- 
nation of experimental reaction rates and high-level elec- 
tronic structure calculations [49j . The mean signed error 
of standard enthalpies of formation is roughly equal and 
opposite to the mean signed error of the corresponding 
atomization energies, whereas the mean absolute errors 
are about the same. The rather complicated methodol- 
ogy of evaluating the standard enthalpies of formation is 



documented in Ref. The only difference between the 
present procedure and the one of Ref. [H is that here we 
employed the fully uncontracted 6-311++G(3df,3pd) ba- 
sis set for calculating the electronic energies rather than 
the standard (contracted) 6-311++G(3df,3pd) basis [78[, 
because the latter is not well-suited for the resolution of 
the identity used for evaluating e ® x ( conv ) anc j £ ™ [57| . 

We do not have yet a fully self-consistent implemen- 
tation of our local hybrid functional and so evaluate all 
hyper-GGA energies using converged TPSS orbitals. We 
speculate, based on our studies of a previous version of 
the hyper-GGA, that the choice of orbitals would have 
a small effect on the results for enthalpies of formation 
and barrier heights, provided that the parameters of the 
functional are optimized for each choice. 

The optimized values of the hyper-GGA parameters 
obtained in this manner are as follows: A — 2.74, B = 
132, C = 0.940, D = 6.13, and E = 8.02. Although 
the G2/97 and BH42/03 training sets differ markedly in 
size, their equal weighting in Eq. (|2"9"j) does not cause a 
significant bias in favor of the data included in either test 
set because both sets are sufficiently representative of the 
types of molecules/reactions they contain. 

In Table HI we compare the performance of our PSTS 
hyper-GGA with that of simpler functionals constructed 
by the method of constraint satisfaction with no or min- 
imal empiricism. The functionals not already introduced 
in Section |T] are: HFx/TPSSc is the Hartree-Fock ex- 
change with TPSS correlation, PBEh is the global hybrid 
PBE functional with a = 0.25 0, H H, and TPSSh 
is the global TPSS hybrid with a = 0.10 [H|. Among 
the approximations included in Table|Tl the PSTS hyper- 
GGA is overall the most accurate (has the smallest Sw)- 
To put these results into perspective, we point out that 
the mean experimental atomization energy for the G2/97 
set of molecules is 478.5 kcal/mol, while the mean ex- 
perimental barrier height for the BH42 set of reaction 
barriers is 14.0 kcal/mol. Note that although the PSTS 
hyper-GGA has been trained on G2/97, it has an even 
better performance for the larger G3/99 test set [44[ of 
223 molecules (including COF2), for which the mean ex- 
perimental atomization energy is ~1180 kcal/mol. 

Table Q] also shows results for two simplified hyper- 
GGAs labeled "hyper-GGA (conv.)" and "hyper-GGA 
(a = ai)" . The former differs from the PSTS hyper-GGA 
in that it uses the conventional exact-exchange energy 
per electron e ° x ( conv ' [ n pl ace f the gauge-transformed 
e° x of Eq. (|12p . The five empirical parameters of this 
form, determined by minimizing Sw, are as follows: 
A = 3.14, B = 146, C = 0.930, D = 5.17, and E = 9.49. 
This form performs slightly worse than the PSTS hyper- 
GGA, which is consistent with our argument that local 
hybrids should combine exact and semi-local exchange in 
the same gauge. The second simplified hyper-GGA uses 
e° x in the TPSS gauge but has ai set to zero, so its mix- 
ing fraction a = a\ depends only on parameters A and 
B (whose optimized values are A = 2.77 and B = 545). 
Because this a is close to 0.1 for nearly all relevant val- 
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TABLE I: Relative performance of the PSTS and two simplified hyper-GGAs for the G2/97 and G3/99 test sets of standard 
enthalpies of formation (148 and 223 molecules, respectively) and the BH42/03 set of 42 reaction barriers. The test sets are 
described in Section [V] G2/97 is a subset of G3/99. ME is the mean error (signed) relative to reference values, MAE is the 
mean absolute error. Sw is defined by Eq. (|29[l . All calculations were carried out with a development version of the GAUSSIAN 
program [8(J using the fully uncontracted 6-311++G(3df,3pd) basis set. All values are in kcal/mol. For comparison, the mean 
atomization energies are 478.5 kcal/mol for G2/97 and 1180 kcal/mol for G3/99, while the mean reaction barrier height is 14.0 
kcal/mol (BH42/03). 1 kcal/mol « 0.0434 eV « 0.00159 hartree. 



G2/97 BH42/03 G3/99 



Method 


ME 


MAE 


ME 


MAE 


Sw 


ME 


MAE 


Max. ( + ) 


Min. (-, 


Hartree-Fock 


148.0 


148.0 


12.5 


12.8 


80.4 


211.0 


211.0 


581.0 (CsHis) 


-0.6 (BeH) 


HFx/TPSSc 


25.2 


27.7 


4.2 


5.8 


16.7 


27.7 


30.3 


130.0 (Os) 


-23.1 (Si 2 H 6 ) 


LSDA" 


-83.4 


83.4 


-18.1 


18.1 


50.7 


-121.2 


121.2 


0.6 (Li 2 ) 


-344.3 (CioH 8 ) 


PBE GGA 


-16.1 


16.9 


-9.7 


9.7 


13.3 


-21.6 


22.1 


10.8 (Si 2 H 6 ) 


-78.9 (CioHg) 


TPSS meta-GGA 


-5.7 


6.4 


-8.4 


8.4 


7.4 


-6.0 


6.5 


15.0 (SiF 4 ) 


-23.9 (C1F 3 ; 


PBEh 


-2.6 


5.0 


-4.7 


4.7 


4.8 


-5.0 


6.8 


20.5 (SiF 4 ) 


-35.8 (CioHg) 


TPSSh 


-1.9 


4.4 


-6.6 


6.6 


5.5 


-1.6 


4.1 


20.8 (SiF 4 ) 


-18.0 (Si 2 H 6 ; 


PSTS hyper-GGA 6 


-0.6 


4.7 


-1.2 


2.1 


3.4 


-0.2 


4.5 


23.2 (SiF 4 ) 


-19.1 (Si 2 H 6 ) 


Hyper-GGA (conv.) 6 


-0.8 


5.6 


-1.2 


2.3 


3.9 


-0.2 


5.5 


28.7 (SiF 4 ) 


-21.5 (Si 2 H 6 ) 


Hyper-GGA (a = ai) 6 


0.2 


4.5 


-6.6 


6.6 


5.5 


1.5 


4.7 


26.5 (SiF 4 ) 


-19.0 (Si 2 H 6 ) 



"Using the Perdew-Wang representation [l5l | of e^ SDA (r s , f ). 
b All hyper-GGA energies were evaluated using self-consistent 
TPSS orbitals. 



ues of it(r), the performance of the ai-only hyper-GGA 
is very similar to that of TPSSh (a global hybrid with 
a = const = 0.1): atomization energies are more accu- 
rate than from the TPSS meta-GGA, but there is little 
improvement for reaction barriers. 



VI. CONCLUSIONS 

We have argued that semilocal density functionals 
work in many cases because of a justified cancellation of 
errors between exchange and correlation that occurs in 
identifiable "normal" regions. This enables one to con- 
struct nonempirical semilocal functionals (such as LSD, 
PBE GGA, and TPSS meta-GGA) on the first three 
rungs of a ladder of increasingly sophisticated approxima- 
tions. The fourth or hyper-GGA rung requires empirical 
parameters to balance the full nonlocality of correlation 
against that of exact exchange, since known exact con- 
straints say nothing about this balance. (At least this 
is so when we consider only exact constraints on the in- 
tegrated exchange-correlation energy; a correlation fac- 
tor [83[, applied to a nearly-exact exchange hole density 
of an inhomogeneous system [84| . might work or might 
not without empiricism). Our PSTS hyper-GGA is a 
local hybrid, mixing in a fraction a of exact exchange lo- 
cally, according to Eqs. 0, HU), and (25]). Because 
our a tends to one in the high-density limit, our hyper- 
GGA has full exact exchange. It is also one-electron self- 
interaction- free and size-consistent. 

The small relative error (typically of order 1% or less) 
of the TPSS meta-GGA for the atomization energies of 
molecules suggests that these quantities are dominated 
by contributions from normal regions. The much larger 



relative error (typically of order 60%) of the TPSS meta- 
GGA for the barrier heights to chemical reactions sug- 
gests that those quantities contain substantial contribu- 
tions from abnormal regions. Thus we fit our hyper-GGA 
parameters simultaneously to atomization energies (pre- 
serving and slightly improving them, mainly via our local 
mixing fraction ai) and to barrier heights (substantially 
improving them, mainly through our local mixing frac- 
tion 02). 

We have invoked five empirical parameters 
(A, B, C, D, E), probably close to the minimum possible 
number since there are four kinds of abnormal regions. 
But, as more empirical parameters are introduced, there 
is a graver danger of "overfitting" any given limited data 
set. A fit that is unjustifiably good can worsen results 
for systems and properties very different from those that 
have been fitted. We have tried to minimize this danger 
by fitting to forms which take into account known exact 
constraints, physical insights, and paradigm examples. 

To construct our PSTS hyper-GGA, we satisfy exact 
constraints on the density functional for the exchange- 
correlation energy. Of all the constraints possible for 
a hyper-GGA, the only one we have not tried to sat- 
isfy is the non-zero limit (69l | for the correlation energy 
under one-dimensional density scaling to the true two- 
dimensional limit. In this limit, our PSTS hyper-GGA 
correlation energy tends to zero (as in the TPSS meta- 
GGA); the correct non-zero limit seems too complicated 
to incorporate in any simple way, and not very rele- 
vant to most physical systems. Moreover, two exact 
constraints that are guaranteed for all possible densities 
at the TPSS meta-GGA level are no longer so guaran- 
teed in local hybrids like PSTS: the non-positivity of the 
correlation energy and the Lieb-Oxford lower bound on 
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the exchange-correlation energy. Of course, the PSTS 
exchange-correlation energy is negative for all possible 
densities. 

Apart from the need for empirical parameters, con- 
straint satisfaction is the same method used earlier to 
construct the PBE GGA 16] and the TPSS meta- 
GGA [18]. In the hyper-GGA case, we satisfy additional 
exact constraints via a careful interpolation between a 
scmilocal exchange energy density in a normal region 
and the exact exchange energy density in an abnormal 
region. Our method combines (and goes beyond) some 
of the best formal features of other methods, including 
the size-consistency of methods that model the exchange- 
correlation hole density @, HE HE HII an d the proper 
concern for the high-density (or weakly-interacting) and 
low-density (or strongly-interacting) limits of methods 
that model the inte gran d W a [n] of the coupling-constant 
integral (H M, MM, M 



E KC [n] = / daW a [n]. 
Jo 



(30) 



Of course, our PSTS hyper-GGA has its own hole model 
[namely, Eq. (j!5p ] and its own coupling-constant inte- 
grand [68, 85] defined by 



W a [n] = -^-(a 2 £ xc [?i 1/Q ]). 



(31) 



The behavior of our abnormality index or local exact- 
exchange mixing fraction a(r), as illustrated in Figs.[3HSl 
is qualitatively reasonable: For an isolated many-electron 
atom of fixed electron number (Fig. [4} it is small (« 0.1) 
in the valence region but larger (« 0.3) in the higher- 
density core, and it gradually approaches 1 in the rapidly- 
varying exponential tail. For a stretched molecule, it is 
atomic-like around an atom of weakly-fluctuating elec- 
tron number (Fig. [5] for spin-restricted stretched F 2 ), but 
much larger (« 0.4) in the relevant valence region of an 
atom of strongly-fluctuating electron number (Fig. [6] for 
stretched Ne^"). A version of the last effect is what raises 
and improves the barrier heights in Table HI 

Figure 0] shows a slow increase of a in the core as the 
atomic number Z increases from 12 (Mg) to 14 (Si 2+ ) to 
20 (Ca 8+ ). At large Z, there is a uniform density scaling 
to the high-density limit in the core, so that a — » 1 as 
Z — » oo. The slowness of the approach to 1 is not a prob- 
lem, since TPSS meta-GGA (unlike LSDA) exchange is 
accurate for localized core-electron densities. 

Figure [5] shows a seemingly unphysical cusp in a at the 
bond center, where for the symmetric stretched molecule 
F2 the reduced density gradient s — > and thus our 
a — > (although the approach to cannot be plotted 
on the scale of this figure). This cusp is a consequence of 
our Eqs. (flT)|) and (flUj) . since u — > 00 as s — *■ 0. By recov- 
ering semilocality in the limit of slowly-varying density, 
we have introduced an artifact in a small volume around 
the bond center. Stretched NeJ in Fig. [5] also displays a 
bond-center cusp. Even the TPSS meta-GGA exchange 



by itself can have a bond-center artifact [72| . So far as 
we know, these artifacts are harmless. 

The ideas in this paper are general ones that may be 
generally useful. On the other hand, our way of imple- 
menting them is far from unique. In future work, we 
hope to continue to refine and test this approach to the 
hyper-GGA including its self-consistent implementation. 

APPENDIX A: EXPRESSION FOR £ G L2TPSS 

The Gorling-Levy second-order limit of the TPSS cor- 
relation energy per electron e GL2TPSS j g gj ven by 



£ GL2TPSS _ £ GL2rcvPKZB 



^ _|_ ^ 6 GL2rcvPKZB ^ T W y 



where d — 2.8 hartree 1 is a constant and 

GL2revPKZB 



(Al) 



OL2PBE 



(iT> n l) V "T> V "i) 



1 



[1 + C(C,0](— ) 2 E-e 2PBE - (A2) 



In Eq. (JA2"1) . £ gl2Pbe ig the Gorling-Levy n m [ t f t h e 
PBE correlation energy per electron. It is obtained by 
replacing Ar with r in the A — ► 00 uniform density scaling 
limit of the PBE correlation energy per electron [l|| : 



GL2PBE 



-jef) 3 In 



XS 2 /4> 2 + (x s2 /0 2 ) 2 



in which 7 = (1 — ln2)/7r 2 



1 



<X0 = 2 (i + C) 2/3 + (i-C) 2/;; 



(A3) 



(A4) 



s = \\7n\/2nkp is the reduced density gradient of 
Eq. (HOD, k F = (Stt 2 ™) 1 / 3 , and X = (P/l)c 2 e- UJ /'> « 
0.72161, where c = {iir 2 /^) 1 / 3 , (5 = 0.066725, and 
lo = 0.046644 111. The spin-dependent function e^ 2PBE 
is defined H,l72| as 



-GL2PBE 
t c,cr 



max[e^ L2PBE (n CT ,0,Vn CT ,0), 



,GL2PBE 



(n h n u Vn T , Vnj.)]. (A5) 



The function C(£,£), where £ is the spin-polarization 
of Eq. © and £ = |VC|/2Jfe F , is given by Eq. (14) of 
Ref. El 



APPENDIX B: EVALUATION OF f(v) 

The exponent 1/(1 — v) — l/(v — C) F appearing in 
Eq. (f2"3")l tends to +c>o as v approaches 1 from below. 
In order to avoid numerical overflow problems with the 
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exponential term of f(v) in the interval C < v < 1, 
Eq. ([2"3"|) may be rewritten as follows. Define 

*=(T3S)F' pc= (^cp- (B1) 

Then, inside the interval C < v < 1, 

!— -— — , if pi <p c 
1 + eP f- P r , (B2) 

1 + e Pc- Pl ■ if ^<Pi 

which is robust because all exponents are non-positive. 
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